﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Drawing;
using System.Drawing.Imaging;

namespace FourierTransform
{
    public struct Complex
    {
        public float RealPart, ImaginePart;

        public static Complex FromReal(float num)
        {
            Complex rs;
            rs.RealPart = num;
            rs.ImaginePart = 0.0f;
            return rs;
        }

        public float GetLength()
        {
            return (float)Math.Sqrt(RealPart * RealPart + ImaginePart * ImaginePart);
        }

        public static Complex operator +(Complex lhs, Complex rhs)
        {
            Complex rs;
            rs.RealPart = lhs.RealPart + rhs.RealPart;
            rs.ImaginePart = lhs.ImaginePart + rhs.ImaginePart;
            return rs;
        }

        public static Complex operator *(Complex lhs, Complex rhs)
        {
            Complex rs;
            rs.RealPart = lhs.RealPart * rhs.RealPart - lhs.ImaginePart * rhs.ImaginePart;
            rs.ImaginePart = lhs.RealPart * rhs.ImaginePart + lhs.ImaginePart * rhs.RealPart;
            return rs;
        }

        public static Complex operator *(Complex lhs, float rhs)
        {
            Complex rs;
            rs.RealPart = lhs.RealPart * rhs;
            rs.ImaginePart = lhs.ImaginePart * rhs;
            return rs;
        }
    }

    public delegate void NotifyEvent();
    public delegate bool BoolEvent();
    public delegate void ProgressChangedEvent(int value, int max);

    public class FourierImage
    {
        public List<Complex> ValuesR, ValuesG, ValuesB;
        public int Width, Height;

        public static FourierImage FromImage(Image img, bool log, BoolEvent shouldCancel, ProgressChangedEvent progChange)
        {
            FourierImage fimg = new FourierImage();
            fimg.Width = img.Width;
            fimg.Height = img.Height;
            Bitmap bit = new Bitmap(img);
            BitmapData data = bit.LockBits(new Rectangle(0, 0, bit.Width, bit.Height), System.Drawing.Imaging.ImageLockMode.ReadOnly,
                System.Drawing.Imaging.PixelFormat.Format32bppRgb);
            List<Complex> tmpR = new List<Complex>(img.Width * img.Height);
            List<Complex> tmpG = new List<Complex>(img.Width * img.Height);
            List<Complex> tmpB = new List<Complex>(img.Width * img.Height);
            unsafe
            {
                for (int i = 0; i < fimg.Height; i++)
                {
                    int* pix = (int*)((int)data.Scan0 + (int)(data.Stride) * i);
                    for (int j = 0; j < fimg.Width; j++)
                    {
                        if (log)
                        {
                            tmpR.Add(Complex.FromReal((float)Math.Log(((pix[j] & 0xFF0000) >> 16) / 255.0f)));
                            tmpG.Add(Complex.FromReal((float)Math.Log(((pix[j] & 0x00FF00) >> 8) / 255.0f)));
                            tmpB.Add(Complex.FromReal((float)Math.Log(((pix[j] & 0x0000FF)) / 255.0f)));
                        }
                        else
                        {
                            tmpR.Add(Complex.FromReal(((pix[j] & 0xFF0000) >> 16) / 255.0f));
                            tmpG.Add(Complex.FromReal(((pix[j] & 0x00FF00) >> 8) / 255.0f));
                            tmpB.Add(Complex.FromReal(((pix[j] & 0x0000FF)) / 255.0f));
                        }
                    }
                }
            }
            bit.UnlockBits(data);
            FourierTransform.FourierTransform2D(tmpR, fimg.Width, fimg.Height, out fimg.ValuesR, shouldCancel,
                (int value, int max) => { progChange(value, max * 3); });
            FourierTransform.FourierTransform2D(tmpG, fimg.Width, fimg.Height, out fimg.ValuesG, shouldCancel,
                (int value, int max) => { progChange(value + max, max * 3); });
            FourierTransform.FourierTransform2D(tmpB, fimg.Width, fimg.Height, out fimg.ValuesB, shouldCancel,
                (int value, int max) => { progChange(value + max * 2, max * 3); });
            return fimg;
        }
        public Bitmap GetSpectrumBitmap(float ample)
        {
            Bitmap bit = new Bitmap(Width, Height);
            BitmapData data = bit.LockBits(new Rectangle(0, 0, bit.Width, bit.Height), ImageLockMode.WriteOnly,
                 PixelFormat.Format32bppRgb);
            unsafe
            {
                for (int i = 0; i < Height; i++)
                {
                    int* pix = (int*)((int)data.Scan0 + (int)(data.Stride) * i);
                    int ni = i - Height / 2;
                    if (ni < 0)
                        ni = ni + Height;
                    for (int j = 0; j < Width; j++)
                    {
                        int nj = j - Width / 2;
                        if (nj < 0)
                            nj = nj + Width;
                        int r = Math.Min((int)Math.Round(Math.Log(ValuesR[ni * Width + nj].GetLength() + 1.0) * 255 * ample), 255);
                        int g = Math.Min((int)Math.Round(Math.Log(ValuesG[ni * Width + nj].GetLength() + 1.0) * 255 * ample), 255);
                        int b = Math.Min((int)Math.Round(Math.Log(ValuesB[ni * Width + nj].GetLength() + 1.0) * 255 * ample), 255);

                        pix[j] = Color.FromArgb(r, g, b).ToArgb();
                    }
                }
            }
            bit.UnlockBits(data);
            return bit;
        }

        public Bitmap ToOriginalBitmap()
        {
            return ToOriginalBitmap(false, null, null);
        }

        public Bitmap ToOriginalBitmap(bool exp, BoolEvent shouldCancel, ProgressChangedEvent progChange)
        {
            Bitmap bit = new Bitmap(Width, Height);

            List<Complex> tmpR;
            List<Complex> tmpG;
            List<Complex> tmpB;

            FourierTransform.InverseFourierTransform2D(ValuesR, Width, Height, out tmpR, shouldCancel,
                (int value, int max) => { if (progChange != null) progChange(value, max * 3); });
            FourierTransform.InverseFourierTransform2D(ValuesG, Width, Height, out tmpG, shouldCancel,
                (int value, int max) => { if (progChange != null) progChange(value + max, max * 3); });
            FourierTransform.InverseFourierTransform2D(ValuesB, Width, Height, out tmpB, shouldCancel,
                (int value, int max) => { if (progChange != null) progChange(value + max * 2, max * 3); });

            BitmapData data = bit.LockBits(new Rectangle(0, 0, bit.Width, bit.Height), ImageLockMode.WriteOnly,
                 PixelFormat.Format32bppRgb);
            unsafe
            {
                for (int i = 0; i < Height; i++)
                {
                    int* pix = (int*)((int)data.Scan0 + (int)(data.Stride) * i);
                    for (int j = 0; j < Width; j++)
                    {
                        if (exp)
                        {
                            pix[j] = Color.FromArgb((int)(Math.Exp((int)(tmpR[i * Width + j].RealPart * 255.0f) & 0xFF)),
                                                    (int)(Math.Exp(((int)(tmpG[i * Width + j].RealPart * 255.0f) & 0xFF))),
                                                    (int)(Math.Exp(((int)(tmpB[i * Width + j].RealPart * 255.0f) & 0xFF)))).ToArgb();
                        }
                        else
                        {
                            pix[j] = Color.FromArgb(((int)(tmpR[i * Width + j].RealPart * 255.0f) & 0xFF),
                                     ((int)(tmpG[i * Width + j].RealPart * 255.0f) & 0xFF),
                                     ((int)(tmpB[i * Width + j].RealPart * 255.0f) & 0xFF)).ToArgb();
                        }
                    }
                }
            }
            bit.UnlockBits(data);
            return bit;
        }
    };

    public static class FourierTransform
    {

        public static void FourierTransform1D(List<Complex> inputValues, List<Complex> outValues)
        {
            float[] realIn = new float[inputValues.Count];
            float[] imagineIn = new float[inputValues.Count];
            for (int i = 0; i < inputValues.Count; i++)
            {
                realIn[i] = inputValues[i].RealPart;
                imagineIn[i] = inputValues[i].ImaginePart;
            }
            float[] realOut = new float[inputValues.Count];
            float[] imagineOut = new float[inputValues.Count];
            unsafe
            {
                Fourier.FourierTransform(inputValues.Count,
                    (float*)System.Runtime.InteropServices.Marshal.UnsafeAddrOfPinnedArrayElement(realIn, 0),
                    (float*)System.Runtime.InteropServices.Marshal.UnsafeAddrOfPinnedArrayElement(imagineIn, 0),
                    (float*)System.Runtime.InteropServices.Marshal.UnsafeAddrOfPinnedArrayElement(realOut, 0),
                    (float*)System.Runtime.InteropServices.Marshal.UnsafeAddrOfPinnedArrayElement(imagineOut, 0));
            }
            outValues.Clear();
            Complex divider = Complex.FromReal(1.0f / inputValues.Count);
            for (int i = 0; i < inputValues.Count; i++)
            {
                Complex c;
                c.RealPart = realOut[i];
                c.ImaginePart = imagineOut[i];
                c *= divider;
                outValues.Add(c);
            }

        }

        public static void InverseFourierTransform1D(List<Complex> inputValues, List<Complex> outputValues)
        {
            float[] realIn = new float[inputValues.Count];
            float[] imagineIn = new float[inputValues.Count];
            for (int i = 0; i < inputValues.Count; i++)
            {
                realIn[i] = inputValues[i].RealPart;
                imagineIn[i] = -inputValues[i].ImaginePart;
            }
            float[] realOut = new float[inputValues.Count];
            float[] imagineOut = new float[inputValues.Count];
            unsafe
            {
                Fourier.FourierTransform(inputValues.Count,
                    (float*)System.Runtime.InteropServices.Marshal.UnsafeAddrOfPinnedArrayElement(realIn, 0),
                    (float*)System.Runtime.InteropServices.Marshal.UnsafeAddrOfPinnedArrayElement(imagineIn, 0),
                    (float*)System.Runtime.InteropServices.Marshal.UnsafeAddrOfPinnedArrayElement(realOut, 0),
                    (float*)System.Runtime.InteropServices.Marshal.UnsafeAddrOfPinnedArrayElement(imagineOut, 0));
            }
            outputValues.Clear();
            Complex divider = Complex.FromReal(1.0f / inputValues.Count);
            for (int i = 0; i < inputValues.Count; i++)
            {
                Complex c;
                c.RealPart = realOut[i];
                c.ImaginePart = -imagineOut[i];
                outputValues.Add(c);
            }
        }

        public static bool FourierTransform2D(List<Complex> values, int width, int height, out List<Complex> result)
        {
            return FourierTransform2D(values, width, height, out result, null, null);
        }

        public static bool FourierTransform2D(List<Complex> values, int width, int height, out List<Complex> result,
            BoolEvent shouldCancel, ProgressChangedEvent progChange)
        {
            List<Complex> tmpX = new List<Complex>(width);
            List<Complex> output = new List<Complex>(width);
            result = new List<Complex>(width * height);
            for (int i = 0; i < height; i++)
            {
                tmpX.Clear();
                for (int j = 0; j < width; j++)
                {
                    tmpX.Add(values[i * width + j]);
                }
                FourierTransform1D(tmpX, output);
                for (int j = 0; j < width; j++)
                {
                    result.Add(output[j]);
                }
                if (shouldCancel != null && shouldCancel())
                    return false;
                if (progChange != null)
                    progChange(i, height + width);
            }
            List<Complex> tmpY = new List<Complex>(height);
            for (int i = 0; i < width; i++)
            {
                tmpY.Clear();
                for (int j = 0; j < height; j++)
                {
                    tmpY.Add(result[j * width + i]);
                }
                FourierTransform1D(tmpY, output);
                for (int j = 0; j < height; j++)
                {
                    result[j * width + i] = output[j];
                }
                if (shouldCancel != null && shouldCancel())
                    return false;
                if (progChange != null)
                    progChange(i + height, height + width);
            }
            return true;
        }

        public static bool InverseFourierTransform2D(List<Complex> input, int width, int height, out List<Complex> result)
        {
            return InverseFourierTransform2D(input, width, height, out result, null, null);
        }

        public static bool InverseFourierTransform2D(List<Complex> input, int width, int height, out List<Complex> result,
            BoolEvent shouldCancel, ProgressChangedEvent progChange)
        {
            result = new List<Complex>(input);
            List<Complex> output = new List<Complex>(height);
            List<Complex> tmpY = new List<Complex>(height);
            for (int i = 0; i < width; i++)
            {
                tmpY.Clear();
                for (int j = 0; j < height; j++)
                {
                    tmpY.Add(input[j * width + i]);
                }
                InverseFourierTransform1D(tmpY, output);

                for (int j = 0; j < height; j++)
                {
                    result[j * width + i] = output[j];
                }
                if (shouldCancel != null && shouldCancel())
                    return false;
                if (progChange != null)
                    progChange(i, width + height);
            }

            List<Complex> tmpX = new List<Complex>(width);

            for (int i = 0; i < height; i++)
            {
                tmpX.Clear();
                for (int j = 0; j < width; j++)
                {
                    tmpX.Add(result[i * width + j]);
                }
                InverseFourierTransform1D(tmpX, output);
                for (int j = 0; j < width; j++)
                {
                    result[i * width + j] = output[j];
                }
                if (shouldCancel != null && shouldCancel())
                    return false;
                if (progChange != null)
                    progChange(i + height, width + height);
            }
            return true;

        }
    }
}
